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An interaction between the vacuum energy and dark matter is an intriguing possibility which 
may offer a way of solving the cosmological constant problem. Adopting a general prescription 
for momentum exchange between the two dark components, we reconstruct ot{a), the temporal 
evolution of the coupling strength between dark matter and vacuum energy, in a nonparametric 
Bayesian approach using combined observational data sets from the cosmic microwave background, 
supernovae and large scale structure. An evolving interaction between the vacuum energy and dark 
matter removes some of the tensions between different data sets. However, it is not preferred over 
ACDM in the Bayesian sense, as improvement in the fit is not sufficient to compensate for the 
increase in the volume of the parameter space. 


PACS numbers: 98.80. -k, 98.80.Es 

I. INTRODUCTION 

The discovery of cosmic acceleration [1, 2] has inspired 
the development of a wide range of dark energy (DE) and 
modified gravity models. The simplest DE candidate, the 
cosmological constant A, is a parameter in General Rela¬ 
tivity and is consistent with all current observations [3]. 
The main problem with A is not why it has a particular 
value, but the fact that the vacuum energy contribution 
to A is sensitive to the ultraviolet cutoff scale and re¬ 
quires a technically unnatural fine-tuning [4], which is 
the long-standing cosmological constant problem (GGP) 
[5] . Most of the dynamical DE and modified gravity mod¬ 
els proposed in the literature do not offer a solution to 
this old problem. The CCP would be surmountable if 
there was a dynamical mechanism by which vacuum en¬ 
ergy could decay from its initially large value and settle 
at an attractor consistent with the observed value of A. 
While a full theory of the quantum vacuum that would 
contain such a mechanism does not exist, this idea has 
motivated phenomenological models of decaying vacuum 
energy [6-12]. 

One way to have a nonconstant vacuum energy is to in¬ 
troduce a new dynamical degree of freedom, e.g., a scalar 
field. An alternative approach, which avoids explicitly in¬ 
troducing new degrees of freedom while still preserving 
the general covariance of the evolution of cosmological 
perturbations, is to allow for the exchange of momen¬ 
tum between the vacuum energy and other species [13]. 
Both the minimally [14-16] and nonminimally [17-19] 
coupled quintessence models can be cast in this general 
framework. In the former case, the vacuum (the poten¬ 
tial energy) exchanges energy with the kinetic energy of 
the scalar field. In the latter, there is an additional ex¬ 
change with matter. In this sense, one can say that a 
time-dependent vacuum energy is necessarily interacting. 
Since additional interactions in the visible matter sector 
are strongly constrained while the nature of dark mat¬ 
ter (DM) is largely unknown, we will consider models in 
which the vacuum interacts only with DM. 


In this paper we adopt the phenomenological model of 
vacuum energy evolution introduced in [13] which avoids 
dealing with explicit additional degrees of freedom. The 
vacuum equation of state is by definition equal to — 1, but 
the vacuum energy density V is allowed to vary because 
of the interaction with DM, = —Qp,- 

II. THE MODEL 

In the interacting vacuum energy model the back¬ 
ground DM and vacuum energy densities obey continuity 
equations 

Pdm + 3i7pdm = -Q-, V = Q , (1) 

where Q denotes the energy transfer between DM and 
vacuum energy. An arbitrary energy transfer Q can in 
principle reproduce an arbitrary background cosmology, 
with energy density p = pdm + V and pressure P = —V 
[13]. This reduces to AGDM when Q = 0. 

To calculate the linear perturbations, we need to spec¬ 
ify a covariant form of the energy-momentum trans¬ 
fer 4-vector. Following [20, 21], we assume that the co¬ 
variant interaction is parallel to the 4-velocity of DM, 
= Qu^dm)- There are other choices [13], but this 
covariant form for the interaction means that there is 
no momentum transfer in the DM comoving-orthogonal 
frame and hence the DM particles follow geodesics, as 
in AGDM. Although the interacting vacuum does allow 
inhomogeneous perturbations, one can always choose a 
frame in which the vacuum is spatially homogeneous, 
SV = 0. For geodesic DM this coincides with comoving- 
orthogonal frame [20, 21] and the perturbation equations 
are particularly simple in this gauge. Note that in this 
case the energy transfer is a potential flow, = — V^U, 

and therefore the matter velocity must be irrota- 

tional [22]. The DM density contrast evolves according 
to 

5dm = -^+ 2^5dm , (2) 

Pdm 
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where the divergence of the matter 3-velocity is given 
by the extrinsic curvature of the metric, "d = h/2 in the 
comoving-synchronous gauge, and h is the scalar mode 
of metric perturbations [23]. 

We take Q to be of a form inspired by the generalized 
Chaplygin gas model [24, 25], Q = 2>aHp^mV/{pdra + V), 
where a is a dimensionless coupling parameter [13]. 
This form naturally reproduces a conventional matter- 
dominated solution at early times, and a vacuum domi¬ 
nated solution at late times, while allowing more general 
evolution in between. Previous studies [20, 21] have con¬ 
strained the interaction assuming a = const. ^ 

In this paper, we make no assumptions about the time 
dependence of a and directly reconstruct it from data us¬ 
ing the nonparametric Bayesian approach introduced in 
[27]. We are thus able to describe a general background 
cosmology and reproduce any equation of state seen in 
quintessence models, but the perturbations are of a re¬ 
stricted form with vanishing sound speed but nonzero 
energy transfer to or from DM, determined by the back¬ 
ground cosmology. Thus our model is degenerate with 
quintessence models in terms of the background cosmol¬ 
ogy, but distinguished by the evolution of perturbations. 

III. THE RECONSTRUCTION METHOD 

We first discretize a(a) into bins ai = a(ai), i = 
1,..., A^, distributed uniformly in the interval [amin, ttmax], 
giving us TV parameters oti that can be fit to data. Any 
representation of a function with a finite number of pa¬ 
rameters is necessarily inaccurate and, in the absence of 
additional priors, the outcome of the fit would directly 
depend on TV. As shown by Crittenden et al. [27], it is 
possible to eliminate the dependence on TV and explicitly 
control the reconstruction bias by adding a prior that 
correlates the nearby bins. 

In the correlated prior approach one assumes that a{a) 
is a smooth function, so that its values at neighboring 
points in a are not entirely independent. More specifi¬ 
cally, a{a) is taken to be a Gaussian random variable 
with a given correlation ^ between its values at a and a', 
^(ja — a'l) = ([Q;(a) — a‘’‘^(a)][Q;(a') — a*’'^(a')]), which is 
nonzero for ja — a'\ below a given “correlation length” 
but vanishes at much larger separations, and is 

a reference fiducial model. Given a particular functional 
form of ^(ja — a'j), one can calculate the corresponding 
covariance matrix for the TV parameters ap. 

/ da/ da'^([a-a'l), (3) 

where A is the bin width. This covariance matrix de¬ 
fines a (Gaussian) prior probability distribution for pa- 


^ Salvatelli et at [26] considered a related dimensionless parameter, 
q = —3opdni/(Pdm -t V), in four redshift intervals. 


rameters ai and its product with the likelihood, ac¬ 
cording to Bayes’ theorem, gives the desired poste¬ 
rior distribution. Thus, the reconstruction amounts to 
finding the minimum of = Xprior + Xdata> where 

Xprior = {ct — C”^ (o — and is some fidu¬ 

cial model. To avoid dependence on we take it to be 
constant and marginalize over its value following [27]. 

We adopt the GPZ form [27, 28] for the correlation 
function, ^(ja — a']) = ^(0)/[l -I- ([a — a'\/acf], where 
ttc determines the correlation length and ^(0) sets the 
strength of the prior. This form has the advantage of 
making it possible to evaluate integrals in Eq. (3) ana¬ 
lytically, while also having a transparent dependence on 
its two parameters. We also note that, as shown in [27], 
the outcome of reconstruction is largely insensitive to the 
particular function chosen to represent ^(ja — a'j). 

The correlation length Qc sets the effective number of 
degrees of freedom allowed by the prior: Ngs ~ (a„iax — 
amin)/ac- As long as TV > Ngg, the reconstructed result 
is independent of TV. Rather than adjusting ^(0) to con¬ 
trol the strength of the prior, we use the variance of the 
mean given by cr? ~ 7r^(0)ac/(aniax — amin) in the limit 
o-c Oniax — Omin [27, 28]. Namely, given Oc, one can 
adjust ^(0) to keep the prior on the uncertainty in the 
mean of a(a) independent of Ngg. In what follows, we 
take (Ta = 0.04 based on the constraint on a constant a 
obtained in [21], and set ag = 0.06. As we show later us¬ 
ing principal component analysis (PGA), this choice of Oc 
effectively separates the signal from noise, thus allowing 
a high-resolution reconstruction with minimal contami¬ 
nation from the noise. 


IV. DATA 

Our data sets include the GMB temperature and po¬ 
larization power spectra from Planck [29] and WMAP9 
[30] respectively; the JLA supernovae sample [31]; the 
BAO measurements of 6dFGRS [32], SDSS DR7 [33], 
BOSS LOWZ [34] and BOSS Lyman-a Forest (LyaF) 
[35]; the redshift space distortion (RSD) measurements 
which probe both the expansion and the growth history, 
namely (TT>y/rs(zd), Pap,/ erg) from BOSS CMASS [36] 
and {A,Fap, fas) from WiggleZ [37], where Fap quanti¬ 
fies the “Alcock-Paezynski” effect [38]. We denote these 
data sets as “RSD (with AP effect),” while we also use 
a “RSD (All)” data set that includes additional RSD 
measurements without the AP effect: 6dFGRS [39], 2dF- 
GRS [40], SDSS LRG [41] and VIPERS [42]. We present 
the constraints from four data combinations as shown in 
Fig. 1. 

The RSD measurements constrain (/erg)^, the product 
of the growth rate and the variance of the total mat¬ 
ter density contrast. Note that the continuity equation 
of DM density contrast in our model [Eq. (2)] is dif¬ 
ferent from that in AGDM and, therefore, the growth 
rate probed by the RSD is no longer simply that of 
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scale factor a 

FIG. 1. Reconstruction of a(a) from four different data combiirations: Data I: CMB+SN+BAO (without LyaF)+RSD (with 
AP effect); Data II: CMB+SN+BAO (without LyaF)+RSD (All); Data III: CMB+SN+BAO (A11)+RSD (with AP effect) and 
Data IV: CMB+SN+BAO (A11)+RSD (All). The best-fit model (central white solid curves) with the 68, 95% CL errors (dark 
and light blue shaded bands) are shown in each panel. We denote the redshift in each panel on the top x-axis. The horizontal 
dashed line denotes the ACDM model. 


the DM component. Namely, in ACDM, the velocity di¬ 
vergence of DM, d = V ■ V, is simply fdmHSdni where 
/dm = dlni5dm/dlna is the DM growth rate. In the in¬ 
teracting vacuum model, Eq. (2) implies that 


t? = - 



Hpdm 


H6dm = -f^HSdm = -UHS^ . 


If Q{ai) ^ 0, the effective DM growth rate fi, the appro¬ 
priately weighted total matter growth rate [43] and 
/dm are different from each other. It can be shown that 


h 


/ 1 Pdm ^ I Ph 

\ fi Pm /b Pm 


(4) 


these cases, the best-fit a{a) model changes sign during 
its evolution, i.e., it is positive at z > 2.1, implying an 
energy transfer from DM to vacuum energy, but negative 
at 0.6 < z < 2.1, implying vacuum decay. At z < 0.6, 
a is consistent with ACDM. The variations at higher 
and lower redshifts compensate to make the deceleration- 
acceleration transition redshift, Zt = 0.6 from Data III, 
close to the ACDM value of Zt = 0.65, which agrees 
with the value extracted from the H{z) data [45] us¬ 
ing a Gaussian prior of Hq = 68 ± 2.8kms“^Mpc“^, 
but is smaller than the value extracted using a prior of 
Hq = 73.8 ± 2.4kms“^Mpc“^ [45]. Again, adding more 
RSD data points (Data IV versus Data III) does not 
change the reconstruction significantly. 


where /b = din 5^/din a is the growth rate of baryons. 
In the absence of an interaction we recover the standard 
result: fs = fi = fh, while, in general, it is (/^cts)^ that 
is measured by RSD. 

V. RESULTS 

We discretize a{a) into 40 bins uniform in a within 
[0.001, 1], and use Monte Carlo Markov chains (MCMC) 
implemented in a modified version of CosmoMC [44] to fit 
them to data along with all other relevant cosmological 
parameters. A total is minimized and the joint pos¬ 
terior probability distribution for all the parameters is 
obtained after the MCMC has converged. 

The best-fit reconstructed models of a{a) (with 68% 
and 95% C.L. errors) are shown in Fig. 1. The ACDM 
fits the observations well when the LyaF BAG measure¬ 
ments are not included (Data I) and the reconstruction 
remains almost unchanged after adding more RSD data 
(Data II). However, when the LyaF data is included, we 
see evidence for an evolving a: Data III (IV) shows a 
1.8(1.9fT) improvement in the fit for the a 0 model. In 


To understand this result, in Fig. 2 we show the the¬ 
oretical predictions for Dy/rs and /^crg for the best-fit 
binned model from Data III rescaled by the correspond¬ 
ing best-fit ACDM model predictions, together with the 
actual measurements. In panel (A), the LyaF measure¬ 
ment (open star), has the smallest error bar and pulls the 
fit down, also making it more consistent with the CMASS 
measurement (filled star). This reduces the BAO by 
4.1. This fit is favored by the fso's as well (panel B), 
further reducing the by 1.7. The remaining data from 
CMB and SN slightly disfavor this model; namely, they 
increase by ^ 2, but not enough to compensate for 
the reduction in x^ from BAO and RSD. 

The improvement in the fit achieved by the binned 
model must be weighed against the increased number of 
degrees of freedom. This can be quantified via the Bayes 
factor, or the ratio of the Bayesian evidences of the in¬ 
teracting model and the ACDM model. For one model to 
be preferred over the other, the Bayes factor should be 
significantly greater than 1. Since the evidence depends 
on the prior assumed for the binned coupling a{ai), we 
follow the method in [47] to examine the dependence of 
the Bayes factor on the the variance in each bin, (Tbin, 
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FIG. 2. Panels (A) and (B): The theoretical prediction of Dy/ra and /^crg by the best-fit a{a) model (solid line), and mea¬ 
surements (data with error bars), rescaled by the best-fit ACDM model. The rescaled BAO measurements are: GdFGRS (filled 
circles), BOSS LOWZ (open circles), SDSS DR7 (filled squares), WiggleZ (open squares), CMASS (filled star), and LyaF (open 
star). The RSD points from WiggleZ (open squares) and GMASS (filled star) are shown in panel (B). The dashed horizontal 
line denotes the ACDM model. Panel (C): The eigenvalues of the covariance matrix obtained using data plus prior for four 
different data combinations, and using prior alone. Panel (D); The first three eigenvectors of the best-fit a{a). See the text for 
more details. 


which controls the strength of the prior. In Fig. 3, we 
plot the logarithm of the evidence ratio (A In E) as well as 
the logarithms of the ratios of the volumes of parameter 
space (AlnV) and the likelihood (AlnL) as a function 
of (Tbin- For CTbin —t 0, the prior effectively forces a(ai) 
to become equal to the fiducial model, which is ACDM, 
and all three ratios approach unity. Increasing crbin allows 
for more freedom in the variation of a(ai) and improves 
the fit. However, as evident from Fig. 3, the improve¬ 
ment in the fit is not sufficient to compensate for the 
increase in the parameter volume. The Bayes factor is 
only marginally greater than I for a limited range of (Tbin 
around the value of O.I, which was the value used in the 
reconstruction in Fig. 1. Thus, we conclude that the in¬ 
teracting model is not preferred over ACDM. 

Even with the lack of strong preference for the evolv¬ 
ing model, it is still interesting to know to what extent 
our reconstruction constrains the evolution of the vac¬ 
uum energy. For instance, one could ask which part of 
the information is informed by data and which part is 
informed by the prior. The errors shown in the recon¬ 
struction in Fig. 1 are highly correlated, making their 
direct interpretation difficult. Principal component anal¬ 
ysis (PCA) is a useful tool that can be used to analyze 
and tune the prior (for applications of PCA to DE stud¬ 
ies, see e.g., [27, 28, 46]). The PCA seeks the orthonor¬ 
mal eigenmodes of the inverse covariance matrix Fa of 
the a bins after marginalizing over other cosmological 
parameters. Namely, Fa = A W, with eigenvectors 
defined by the decomposition matrix W and the eigen¬ 
values given by the elements of the diagonal matrix A. 
Eigenmodes define independent linear combinations of 
the original parameters (a bins) that have uncorrelated 
errors and, thus, are easier to interpret. We can also use 
PCA to compare the eigenmodes with and without the 
prior which, as we explain below, allows us to estimate 
the effective number of degrees of freedom of a(a) con- 



FIG. 3. The logarithms of the ratios of volumes of parameter 
space (AlnV), of likelihoods (AlnL), and of the evidences 
(AlnE) for the binned model and the ACDM as a function 
of the strength of the prior set by the variance in each bin, 

^bin* 


strained by data. Performing the PCA is straightforward, 
since Fa is one of the products of our MCMC calculation. 

The eigenvalues of both the prior and the data-|-prior 
covariance are shown in panel (C) of Fig. 2. There are 
three data eigenmodes that are not affected by the prior, 
le., these three best constrained modes pass the prior 
with almost no penalty. The shapes of these modes are 
shown in panel (D)and we find that they are similar for all 
data combinations. The remaining modes, on the other 
hand, are dominated by the prior. Thus, even though 
our model has many bins of a as parameters, there are 
effectively only three additional degrees of freedom. As 
mentioned earlier, we find that the total is reduced 
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Data I 

Data II 

Data III 

Data IV 


-0.17 ±0.30 

-0.19 ±0.30 

-0.34 ±0.30 

-0.40 ±0.29 

Pi 

0.34 ±0.61 

0.34 ±0.63 

0.47 ±0.65 

0.48 ± 0.64 


-0.25 ± 1.35 

-0.21 ± 1.35 

-0.19 ± 1.34 

-0.10 ± 1.31 


TABLE I. The coefficients (best-fit and 68% C.L. uncertainty) 
of the first three best-determined modes, Pi, of the best-fit 
modeis using different data combinations. 


by ^ 4, which is somewhat greater than expected for a 
model with three additional parameters. 

The eigenvectors provide a natural basis onto which an 
arbitrary a{a) can be expanded, i.e., a(a) = 

Given any a(a), the coefficients /3 can then be found using 
the orthogonality of the modes. The P's corresponding 
to the three best constrained eigenmodes of the best- 
fit models from the four data combinations are listed in 
Table I. Since the uncertainty in the third eigenmode is 
large, adding more modes by relaxing the prior [either by 
reducing ^(0) or ad would not notably change the ht. 


VI. CONCLUSION AND DISCUSSIONS 

We performed a high-resolution reconstruction of a, 
the coupling between DM and vacuum energy, as a func¬ 
tion of the scale factor using the latest observations in¬ 
cluding CMB, SN, BAG and RSD. Our model is de¬ 
generate with standard quintessence models in terms of 
the background cosmology, but is distinguished by the 
growth of perturbations. We found that, when the BAO 
measurement using the BOSS LyaF sample [35] is used. 


an evolving a is mildly favored by the joint data set. In¬ 
terestingly, the best-fit a{a) model changes sign during 
its evolution: a > 0 at higher redshifts, implying an en¬ 
ergy transfer from DM to vacuum energy, while a < 0 
at lower redshifts, corresponding to a decaying vacuum 
energy. A PCA study of our result shows that we have 
extracted three informative eigenmodes from the data. 

The LyaF BAO measurement, which is the best exist¬ 
ing BAO measurement at such a high redshift, is in ten¬ 
sion with the ACDM model at the 2 —2.5(T level. As noted 
in [35, 48], it can be interpreted as favoring a DE compo¬ 
nent with a negative energy density at z ~ 2.3. The RSD 
measurements from BOSS and WiggleZ are also in ten¬ 
sion with ACDM: the RSD measurements favor a lower 
growth rate than the ACDM prediction at low redshifts. 
The interacting vacuum model provides another physical 
interpretation of these tensions if a is allowed to change 
sign during its evolution. We note that extracting BAO 
from LyaF data is a relatively new field, and the cur¬ 
rent measurement could be subject to systematic issues 
[35]. Our method opens a new window into investigation 
of the interacting vacuum model that can be applied to 
improved future datasets as they become available. 
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